A high-order integral solver for scalar problems of diffraction by 
screens and apertures in three dimensional space 

Oscar P. Bruno and Stephane K. Lintner 
^ August 28, 2012 

o 

Abstract 

We present a novel methodology for the numerical solution of problems of diffraction by 
infinitely thin screens in three dimensional space. Our approach relies on new integral formula- 
tions as well as associated high-order quadrature rules. The new integral formulations involve 
weighted versions of the classical integral operators associated with the thin-screen Dirichlet 
and Neumann problems as well as a generalization to the open surface problem of the classical 
Calderon formulae. The high-order quadrature rules we introduce for these operators, in turn, 
resolve the multiple Green function and edge singularities (which occur at arbitrarily close dis- 
tances from each other, and which include weakly singular as well as hyper singular kernels) and 
thus give rise to super-algebraically fast convergence as the discretization sizes are increased. 
When used in conjunction with Krylov-subspace linear algebra solvers such as GMRES, the 
resulting solvers produce results of high accuracy in small numbers of iterations for low and 

i " i high frequencies alike. We demonstrate our methodology with a variety of numerical results for 

screen and aperture problems at high frequencies — including simulation of classical experiments 
such as the diffraction by a circular disc (including observation of the famous Poisson spot), 
interference fringes resulting from diffraction across two nearby circular apertures, as well as 

^ more complex geometries consisting of multiple scatterers and cavities. 

IT) 

1 Introduction 

O 

{N) Diffraction problems involving infinitely thin screens play central roles in the field of wave prop- 

agation: as a noted early example we mention the experimental observation of a bright area in 
J> the shadow of the disc (the famous Poisson spot), which provided one of the earliest confirmations 

of the wave-theory models of light [5, p. xxvhi]. Certainly, problems of diffraction by screens 
continue to impact significantly on a varied range of present day technologies, including wireless 
communications, electronics and photonics, as well as sonic metamaterials, sound transmission, 
non-invasive evaluation and environmental acoustics. Like other wave scattering problems, screen 
problems (mathematically described by open surface problems [25]) can be treated by means of 
numerical methods that rely on approximation of the equations of electromagnetism and acoustics 
over volumetric domains (on the basis of, e.g., finite-difference or finite-element methods) as well as 
methods based on boundary integral equations. Boundary-integral methods do not suffer from the 
well known pollution/dispersion errors [2,14], they require discretization of domains of lower di- 
mensionality than those involved in volumetric methods, and, in spite of the fact that they give rise 
to full matrices, they can be treated efficiently, even for high-frequencies, by means of accelerated 
iterative scattering solvers [4,8,9,23]. Unfortunately, integral methods for open surfaces do not 
give rise, at least in their classical formulations, to Fredholm integral operators of the second-kind, 
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and they suffer, like their volumetric counterparts, from singularities in the vicinity of open edges. 
As a result of these characteristics, both volumetric and boundary-integral numerical methods for 
open surface problems have typically proven computationally expensive and inaccurate. 

Generalizing the two-dimensional open-arc method introduced in [10, 15], this paper presents a 
novel approach for the treatment of open surface diffraction problems in three-dimensional space, 
the benefits of which are two-fold. On one hand, the new method enjoys high-order accuracy: by 
considering weighted versions S w and N w of the single-layer operator S and hypersingular operator 
N, and thus explicitly extracting the solutions' edge singularity, we introduce high-order quadra- 
ture rules (based on the partition of unity method [9] and a combination of polar and quadratic 
changes of variables) which acurately resolve the multiple Green function and edge singularities 
that occur at arbitrarily close distances from each other, and which include weakly singular as 
well as hypersingular kernels. On the other hand, the method gives rise to well-behaved linear 
algebra: as shown in this paper, the composite operator N W S W (which in the two-dimensional case 
was rigorously proven [10, 15] to be a second-kind Fredholm operator) requires very small number 
of iterations when used in conjunction with the linear iterative solver GMRES. In particular, the 
computational times required by our non-accelerated open surface solvers are comparable to those 
required by the corresponding non-accelerated version of the closed surface solver presented in [9] . 
(An extension of the acceleration method [9] to the present context, which is not pursued here, 
does not present difficulties.) Thus, the present methodologies enable solution of the classical open 
surface scattering problems with an efficiency and accuracy comparable to that available for the 
smooth closed-surface counterparts. 

The difficulties associated with boundary-integral methods for open surface problems are of 
course well-known, and significant efforts have been devoted to their treatment. The contributions 
[11, 21] sought to generalize the Calderon relations in the open-surface context as a means to derive 
second-kind Fredholm equations for these problems. In [21] it was shown that the combination NS 
can be expressed in the form I + Tk, where the kernel K(x,y) of the operator Tk has a polar 

singularity of the type O (^ \ x - y \ j ■ This result is not uniform throughout the surface, and it does not 
take into account the singular edge behavior: the resulting operator Tk is not compact (in fact it 
gives rise to strong singularities at the surface edge [15]), and the operator I + T^ is therefore not a 
second-kind operator in any meaningful functional space. When used in conjunction with boundary 
elements that vanish at edges, however, the combination NS can give rise to reduction of iteration 
numbers, as demonstrated in reference [11] through numerical examples concerning low- frequency 
problems. The contribution [11] does not include details on accuracy, and it does not utilize integral 
weights to resolve the solution's edge singularity. A related but different method was introduced 
in [1] which, once again, exhibits small iteration numbers at low frequencies, but which does not 
resolve the singular edge behavior and for which no accuracy studies have been presented. An 
effective approach for regularization of the singular edge behavior in the two-dimensional open-arc 
problem is based on use of a cosine change of variables; see [10] and references therein. In the three 
dimensional case under consideration in this paper, high-order integration rules for the single-layer 
and hypersingular operators were introduced in [25], but these methods have only been applied to 
problems of low frequency, and they have not been used in conjunction with iterative solvers. 

The remainder of this paper is organized as follows. Section 2 defines the Dirichlet and Neumann 
problem on an open surface and it briefly discusses a number of difficulties inherent in classical 
formulations; Section 3 introduces the new weighted operators; and Section 4 provides an outline 
of the Nystrom-based numerical framework on which the solvers are based. The next five sections 
describe the construction of the high-order numerical approximations we use for weighted operators: 
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Sections 5 through 7 decompose the operators into six canonical integral types, while Sections 8 
and 9 provide high-order integration rules for each one of the canonical operators. The selection of 
certain parameters required by our solvers are detailed in Section 10. Finally, numerical results are 
presented in Section 11 which demonstrate the properties of the integral formulations and solvers 
introduced in this paper across a range of frequencies and geometries — including simulation of 
classical experiments such as the diffraction by a circular disc (including observation of the famous 
Poisson spot), interference fringes resulting from diffraction across two nearby circular apertures, 
as well as more complex geometries consisting of multiple scatterers and cavities. 



2 Open-Surface Acoustic Diffraction Problems 

Throughout this paper T denotes a smooth open surface (also called a screen [25]) with a smooth 
edge dT in three-dimensional space. 

2.1 Classical integral equations 

We consider the sound-soft and sound-hard problems of acoustic scattering by the open screen T, 
that is, the Dirichlet and Neumann boundary value problems 

f An + k 2 u = outside T, u\r = f, (sound-soft), . . 

\ Av + k 2 v = outside T, f^| r = g, (sound-hard), ^ ' 

for the Helmholtz equation, where u and v are radiating functions at infinity. 

As is well known, both boundary-value problems are uniquely solvable in adequate functional 
spaces [25]. For r outside T, the solution of the Dirichlet and Neumann problems can be expressed 
as a single-layer potential 

«(r) = ^G fe (r,r>(r')dS' (2) 

and a double-layer potential 

»(r) = / r ^£V)<(S', (3) 
respectively, where Gk denotes the free-space Green's function 

ifc|r-r'| 

Letting S and N denote the classical single-layer and hypersingular operators 

S[u](r) = j G k (ry)fi(r)dS', r on T (5) 

(6) 



and 



N i/ r =hm— / y ' , -^vlr^dS', r on T 

z^o dn r J r dn' r 

the densities \x and v are the unique solutions of the integral equations 

S[/x] = / and (7) 
NH = g. (8) 



3 



The integral operators in (7) and (8) have eigenvalues which accumulate at zero and infinity re- 
spectively, and, thus, solutions of (7) and (8) by means of Krylov-subspace iterative solvers such as 
GMRES generally require large number of iterations. Furthermore, as discussed in Section 2.3, the 
solutions fi and v are singular at the edge of T and thus give rise to low order convergence unless 
such singularities are appropriately taken into account. 



2.2 Calderon formulation in the case of closed surfaces and shortcomings in a 
direct extension to open surfaces 

In the case where the surface T c under consideration is closed (that is, it equals the boundary 
of a bounded set in space), second kind Fredholm equations can be derived either by making 
use of the classical jump relations across the surface of the double-layer potential or the normal 
derivative of the single-layer potential [13], or, alternatively, by relying on the Calderon formula 
which establishes that the combination N C S C of the closed-surface hypersingular operator N c and 
single-layer operator S c can be expressed in the form 

N c S c = -^ + K c , (9) 

where K c is a compact operator in a suitable function space [20]. 

In the case of an open surface, the requirement (1) that the same limit be achieved on both sides 
of the surface prevents the use of discontinuous potentials. And, use the Calderon formula (9) does 
not give rise to a Fredholm equation in the function spaces associated with open-screen problems: 
for example, the composition of N and S is not even defined in the functional framework set forth 
in [25] — since, as shown in [10, 15], the image of the operator S (the Sobolev space H^(T)) is larger 
than the domain of definition of N (the Sobolev space H^(T)). It is interesting to note, further, 
that, as shown in [15] in the two-dimensional case, the image of a constant function has a strong 
edge singularity, 

NS[l](r) = 0(JL), (10) 

where d denotes the distance to the edge — which demonstrates the degenerate character of the 
composite operator NS. 



2.3 Singular edge-behavior 

Although issues related to the singularity of the solutions [i and v of equations (7) and (8) were 
controversial at times [3, 6], the singular character of these solutions is now well documented [7, 16- 
18, 25]. In particular, in reference [16] it was shown that fi and v can be expressed in the forms 

M~^=) z/~Vv / d, (11) 

where ip and ip are infinitely differentiable functions throughout T, up to and including the edge. 
Thus the singular behavior of the solutions to (7) and (8) is fully characterized by the factors d 1 / 2 
and d~ 1 / 2 in equation (11). 



4 



3 Weighted operators and regularized formulation 



In view of the regularity results (11) we introduce a weight u(r) which is smooth, positive and non 
vanishing across the interior of the surface, and which has square-root asymptotic edge behavior: 

u ~ al 1/2 . (12) 

We then define the weighted operators 

S W [ ¥ »] = S[^1, N w [f]=NM, (13) 

-OJ - 

so that for functions / and g that are smooth on T, up to and including the edge dF, the solutions 
of the equations 

SM = f, (14) 

and 

N u ty]=g (15) 

are also smooth throughout the surface. In view of the closed-surface Calderon formula (9), we 
consider the combined operator N W S W and the corresponding equations 

N W S W [¥>] = N W [/] and (16) 

N w S w M=0- (17) 
Note that the solutions of equations (14) and (15) are related to those of (16) and (17): 

»=-, v = u-SM. (18) 

OJ 

As shown in [10, 15] in the case of an open-arc in two dimensions, the combination N^S^ gives 
rise to second-kind integral equations. In particular, the numerical results presented in [10] show 
that equations (16) and (17) require significantly smaller numbers of GMRES iterations than 
equations (14) and (15) to achieve a given residual tolerance. As demonstrated in this paper 
through a variety of numerical examples, a similar reduction in iteration numbers results for three- 
dimensional problems as well. 



4 Outline of the proposed Nystrom solver 
4.1 Basic algorithmic structure 

In order to obtain numerical solutions of the surface integral equations (14)-(17) we introduce an 
open-surface version of the closed-surface Nystrom solver put forth in [9]. This algorithm relies on 

1. A discrete set of nodes N = {ri,i = 1, . . . ,N} on the surface T, which are used for both 
integration and collocation; 

2. High-order integration rules which, using a given discrete set of accurate approximate values 
(<Pi) (resp. (ipi)) of a smooth surface density ip, (pi ~ ^( r «) (resp. ip, ipi ~ ipfa)), pro- 
duce accurate approximations of the quantities S w [y](rj) (resp. N w [V>](rj)), see Sections 4.2 
through 4.4; and 
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Figure 1: Patches, partition of unity and discretization for a disc. Left: the disc is covered by 
an interior patch and two edge patches. Right: partition of unity functions Wf supported on the 
patches. Notice the quadratic refinement along the edges. 



3. The iterative linear algebra solver GMRES [24], for solution of the discrete versions of equa- 
tions (14)-(17) induced by the approximations mentioned in point 2. 

The fact that the same set of Nystrom nodes is used for integration and collocation facilitates 
evaluation of the composite operator N U S W through simple subsequent application of the discrete 
versions of the operators and N w . 

4.2 Partition of unity and Nystrom nodes 

The integration rules mentioned in point 2 in Section 4.1 rely on a decomposition of a given open 
surface (screen) T as a union 



of overlapping patches, including interior patches Vf, q = 1,...,Q%, and edge patches V%, q = 
1, . . . , Qi- For each q, the interior patch V\ (resp. edge patch Pf ) is assumed to be parametrized by 
an invertible smooth mapping = rf(u, v), : Tif — > V\ (resp. x:\ = v\{u, v), n\ : ~H\ — > "Pf) 
defined over an open domain U\ C M 2 (resp. U\ C R 2 f]{v > 0}). Note that, for the q-th edge 
patch, the restriction of the mapping to the set %\ n {v = 0} (which we assume is non-empty for 
q = 1, . . . , Q2) provides a parametrization of a portion of the edge of T; see Figure (1) . Following [9] , 
further, we introduce a partition of unity (POU) subordinated to the set of overlapping patches 
mentioned above. In detail, the POU we use is a set of non-negative functions and defined 
on r, q = 1, . . . , Qi and q = 1, . . . , Q2, such that, for all q, Wf vanishes outside V\ , W$ vanishes 
outside 7^*2 , and the relation 




(19) 




q=l 



q=l 
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holds throughout T. The POU can be used to decompose the integral of any function over the 
surface into a patch-wise sum of the form 

2 Qj 

f f{v')dS' = J2J^ [ j(r q (u,v)) W q (r q (u,v)) J q (u,v)dudv, (20) 

JT j=l q=l JU ) 

where J q (u, v) denote the Jacobian associated with the parametrization r q . At this stage we define 

the set of Nystrom nodes: introducing, for each q, a tensor-product mesh {(u q ' 3 , Vm)} within H q 

(for j = 1, 2), we obtain points = r q (uf 3 , v m ) on the surface V. For every j = 1, 2 and every 

q = l,...,Qj, the set of nodes r\ ,3 m for which the POU function W q associated with the patch V q 
does not vanish 

j^J = { r « : wjftjj > o} (21) 

defines the set of Nystrom nodes on the the patch V q . The overall set N of Nystrom nodes on the 
surface V mentioned in Section 4.1, is given by 

2 Qj 

M=\J\jM q ' j . (22) 

j=l q=l 

Remark 1 . For later reference we introduce the classes of functions 

2?? = {<P\ G C°° (W?) : supp(^) «g W?} , q = 1, . . . , Q 1 , (23) 

V\ = {$\ G C°° (Hi) : supp(^) «g %l} , g = 1, . . . , Q 2 , (24) 

where C°° (Hf) denotes the set of infinitely differentiable functions defined on the open set "H^, 
C°° (%2) denotes the set of functions defined on the set H\ that are infinitely smooth on 'H\ up 
to and including the edge Ji\ fl {v = 0}, and where, for sets A and BCR" the notation A B 
indicates that the closure of A in W 1 is a compact subset of W 1 that is contained in B. 

4.3 Canonical decomposition and high-order quadrature rules 

The high-order numerical quadratures required by step 2. in Section 4.1 are obtained by applying 
the patch-wise decomposition (20) to the weighted integral operators and N w and reducing each 
one of them to a sum of patch-integrals which, as shown in Sections 5 through 7, can be classified 
into six distinct canonical types. For each of these canonical types we construct, in Sections 8 
and 9, spectrally convergent quadrature rules on the basis of the Nystrom points N qj . Suggestions 
concerning selection of sizes of the edge and interior patches, which are dictated on the basis of 
efficiency and accuracy considerations, are put forth in Section 10. 

4.4 Computational implementation and efficiency 

The high-order methods presented in the following sections enable accurate and fast computation 
of each one of the six canonical integral types mentioned in Section 4.3. For added efficiency, 
however, our solver exploits common elements that exist between the various canonical integration 
algorithms, to avoid re-computation of quantities such as the trigonometric functions associated 
with the Green's function, partition of unity functions, integral weights, etc. Additional efficiency 
could be gained by incorporating an acceleration method (see e.g. [9] and references therein) and 
code parallelization. 
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Figure 2: Diffraction by an infinitely thin disc: solution to the Neumann problem for a disc of 
diameter 24A under normal incidence. The famous Poisson spot is clearly visible at the center of 
the shadow area; see also Figure 12. The coloring on the disc represents the values of the surface 
unknown ip. 

5 Canonical singular-integral decomposition of the operator S w 

In view of equation (20), we express the weighted single-layer operator S u , 

S„M(r)= / G k (v,r')^-dS', (25) 
Jv w l r ) 

in the form 

Qi Q 2 

S„ = I>i+E^ ( 26 ) 

g=l g=l 

where 

S q M(v)= Gk M(u,v),r)^ {-J](u,v)WUr q j (u,v))dudv, j = 1,2. (27) 

Jh] w(r?(«,«)) 

In Sections 5.1 and 5.2 the integrals (27) are expressed in terms of canonical integrals of various 
types. 



5.1 Interior patch decomposition 

For an interior patch V\ and for a point r G r \ V\ , the integrand in (27) is smooth and compactly 
supported within the domain of integration %\ — since the weight o>(r) is smooth and nonzero away 
from the edge, and since the POU function W q vanishes outside H\ — and, thus, the integral (27) 
gives rise to our first canonical integral type: 
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Figure 3: Diffraction by an infinitely thin disc: solution to the Dirichlet problem for a disc of 
diameter 24A under incidence parallel to the disc. The coloring on the disc represents the values of 
the surface unknown <j). 



Canonical Integral of Type I 

lf' re9 [cpi] = / (f>i(u,v)dudv , 4>\ G V\ , 
Jul 



(28) 



see Remark 1. For a point r E V\, r = r'(uo,«o) for some (uo,^o) £ H\, on the other hand, the 
integrand of (27) with j = 1 has an integrable singularity at the point (uq,vq) (cf. equation (4)). 
Following [9] we express the kernel as a sum of a localized singular part, and a smooth remainder, 
G k = G« n9 + CT k ea , where 

G ™ 9 = rirG re, = (1 _ Vr)G r k e + jG *n (29) 

Here, G£ e (r, r') = COS |^j r r ,| r ^ and G^ m (r, r') = sm |^ r ,| r ^ denote the real and imaginary parts of the 
kernel G&(r, r'), respectively, and ry r is a smooth function which vanishes outside a neighborhood of 
the point r. As in the previous reference, the collection of all pairs (rj r , 1 — r/ r ) for r E T is called a 
floating partition of unity. The integral that arises as is replaced in (27) by G r ^ g , has a smooth 
integrand which is compactly supported within 7i\; clearly, this is an integral of canonical type 
I. The integral obtained by substituting Gk by G s ™ 9 , on the other hand, gives rise to our second 
canonical type 

Canonical Integral of Type II 

ir n9 l<Pi](u ,v )= [ ^fidudv, ^GVf (30) 
Jh\ l R l 

where for the sake of conciseness, we have set R = r^uo, ^o) — r i( n i v ) and where the point {uq, vq) 
belongs to H\. 
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5.2 Edge-patch decomposition 

The edge singularity on an edge patch V\ is characterized in terms of the asymptotic form (12). In 
what follows we assume, as we may, that on each edge patch,the weight oj is given by an expression 
of the form 

uj(r q 2 (u, v)) = lu%(u, v)y/v, (31) 

where the function uj^u, v) is smooth up to the edge and it does not vanish anywhere along the 
edge. It follows that for an edge patch V\ and for a point r € T \ P 2 ? , the operator <S| defined in 
equation (27) for j = 2 takes the form of an integral of our third canonical type: 



Canonical Integral of Type III 

rwK«»)=/«v)^, fc€i>«. (32) 

Jul V v 

Finally, for an edge patch and for a point r £ V\ we once again use the floating partition of 
unity to decompose the Green function as a sum of a singular and a regular term. The regular term 
results in a canonical integral of Type III, and the singular term gives rise to our fourth canonical 
type: 



n 2 



R 

- '-2 

where R = r^tio^o) — v 2( u -> v )- 



Canonical Integral of Type IV 

ir n9 lM= I Mu \ v)d ^, frev* (33) 



6 Canonical decomposition of the operator N w 

In view of equations (6) and (13) the operator N w is given by the point-wise limit 

N w M(r) ee Jim A jf dGk{r ^ Znr,) ^r')u;(r')dS', r G T. (34) 

Following the open-arc derivation [12, 15, 19] we obtain an adequate expression for this "hypersin- 
gular operator" by taking advantage of the following lemma. 

Lemma 1. The operator N w can be expressed in the form 

N w = Nf, + (35) 

where, denoting the surface gradient with respect to r' by Vp and letting [ • , • ] denote the vector 
product, the operators N£>, NS" and are given by 

r»](r) = fc 2 jT G fc (r,r , )^(r , )a;(r , )n r /.nrd5', (36) 

Tu#](r') = u, 2 (r')VMV>](r') + ^V?V](r'), and (37) 
C[T](r)=p... / [V r G fe (r,r'), [rv,T(r')] ] ■ Ur-^. (38) 
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Proof. See Appendix A. □ 

As shown in Sections 6.1 through 6.3, we may evaluate N w by relying on Lemma 1 and using 
quadratures for various types of canonical integrals. 

6.1 Canonical decomposition of the operator (equation (36)) 
Calling ^2 = ipu 2 we re-express (36) in the form 



N & M(r) = ^/ Gfc (ry) (M^) ^. 
Jr w l r ) 



(39) 



Since u 2 (r) is a smooth function of r throughout T, a construction similar to the one used for (25) 
yields a decomposition of the operator in terms of canonical integrals of types I- IV; see Section 5. 

6.2 Canonical decomposition of the operator %j (equation (37)) 

Making use once again of the POU introduced in Section 4.2, we obtain the decomposition 

Qi Q2 
V? M (r) = £ V? mi) (r) + £ V? Wfl (r ) , 

of the surface gradient. The evaluation of the q-th term in each one of these sums requires the 
calculation of partial derivatives of the form 

d(/>i(u,v) d(f)i{u,v) ^ 



du dv 
^u,v) ^ dMu,v) 



du dv 

for functions <p\ € T>\ and (f>2 € T>\. These partial derivatives can be evaluated efficiently and with 
high-order accuracy by means of the differentiation methods introduced in Sections 8 and 9 below. 
In view of equation (37), use of such high-order rules enables high-order evaluation of the operator 

6.3 Canonical decomposition of the operator (equation (38)) 

It is easy to check that, for a smooth field T, N£T[T] can be evaluated as a linear combination of 
functions of the form D^[^], where is a vector quantity that varies with r but is independent 
of r', where the operator D^f is defined by 

D]fM(r) = P-v. / {V r G fe (r,r') • V} (42) 

and where 4>e are smooth functions. Applying the decomposition (20) to the operator defined in 
equation (42) yields 

Qi Q2 

q=l q=l 
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where 

Vj' q Mr) = p.v.J {v P G fc (r,r?(u, V )) • V} ^0^Wj(r^u,v))J](u,v)dudv. (44) 

We can express the operator D^' 9 as the sum of a hypersingular operator and a weakly singular 
operator whose respective kernels are defined by the split 

V r G fe (r,r')-V = Gf + G%°, Gf = ~ ~ { '^'^ » ( 45 ) 

where the residual kernel G^ s equals a sum of functions which are either smooth or weakly singular 
with singularity [~pi ■ Using the partition-of- unity split embodied in equation (20), the operator 
with kernel G^ s can be expressed in terms of integrals of canonical types I- IV. The hypersingular 
operator with kernel Gj™ on the other hand gives rise to our fifth and sixth canonical types: 

Canonical Integral of Type V 

^i ,pV [<Pi](uo, v ) = p.v. / L Y (f) 1 (u,v)dudv. ^ 

Ju\ l R r 



v f R • V 

^2' PV [(f>2](uo,V ) = p.V. / -^T?-(t>2(u,V 

Jul 



Canonical Integral of Type VI 
R-V dudv (47) 



v 7 ^ 



7 Canonical decomposition of the composite operator N^S 



While the action of the composite operator N W S W on a function <p can be evaluated by producing 
first t/j = S u [(/>] and then evaluating [?/>], both of which can be obtained by the methods described 
in the previous sections, we have found it advantageous in practice to proceed differently, on the 
basis of the expression (35); see Remark 5 for more details. Using the decomposition (35), we first 
evaluate the term N£,S W [<?!>] by means of a direct composition: we compute S u [<p] and then apply 
N£ to the result using the decompositions put forth in Sections 5 and 6.1 respectively. To evaluate 
the second term N^S W [(f)], on the other hand, we first evaluate the quantity 7^S W [0] by expressing 
the surface gradient of S w [</»] required by equation (37) as 

V^M(r)=p., / V*G fe (r,r')^W = £ D^Kr), (48) 

where n(r) and T2(r) denote two orthogonal tangent vectors to T at the point r which vary smoothly 
with r; the corresponding surface gradient of oj 2 can be obtained by direct differentiation of a closed 
form expression, if available, or by means of the differentiation methods put forth in this paper. 
The terms in the sum on the right-hand side of (48) are of the form given in equation (42) with 
V = Tg(r) (£=1,2), and thus can be expressed in terms of the canonical integrals of type I- VI, as 
outlined in Section 6.3. 
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Figure 4: Neumann problem on a spherical cavity of diameter 18A. The coloring on the spherical 
wall represents the values of the surface unknown ip. 

8 High-order evaluation of interior-patch operators 

In this section we describe our algorithms for evaluation of the interior-patch operators introduced 
in the previous sections, namely the integral operators of type I, II and V and the differentiation 
operator (40). To do this, and in accordance with Section 4.2, we assume that the nodes (mJ' 1 ,^ 1 ) 
(uf 1 = Uq 1 + ihu 1 , «m = Vq 1 + mill' 1 , £ = 1, . . . L q , m = 1 . . . Mg), discretize a rectangle that 
contains T~L\. 

8.1 Type I Integral (Regular) 

We evaluate the canonical Type I integral defined in equation (28) by means of a simple trapezoidal 
sum over the grid: as noted in [9] the periodicity (0i is compactly supported) and smoothness of 
the integrand gives rise to super-algebraic convergence in this case. 

8.2 Partial Derivatives 

In view of the smoothness and periodicity of the function (pi, a standard two-dimensional FFT- 
based interpolation scheme based on the evenly spaced grid values 4>^ m yields spectrally convergent 
approximations of the function (f>i(u, v) and its derivatives; our algorithm thus evaluates the deriva- 
tives required in equation (40) by performing a direct term by term differentiation of the resulting 
Fourier representation. 
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Figure 5: Total Field inside a spherical cavity of diameter 18A, Dirichlet Problem. 



8.3 Type II Integral (Singular) 

In order to resolve the singular integrand in equation (30) we utilize the polar change of variables 
introduced in [9]. Defining u(p,9) = uq + pcos9 and v(p,6) = v o + psinO, we obtain 

l q 1 ' sm9 [<f>i](uo,v )= r I q pl [</>x}(u ,v ,9)d9 (49) 
J o 

with 

I^K^fl) = j°° f^pM^dp, (50) 

where 

R=\r q (u ,v )-r q (u(p,9),v(p,9))\, (51) 

and where 

4>i(p, 9) = 4>i(uq + pcos 9, vo + psin^) (52) 

is a smooth function of p and 9 which vanishes for sufficiently large values of p. Since, as noted 
in [9], the ratio ^ is a smooth function of p, the integral I q pl \4>i}{9,uo,vo) defined in (50) can be 
computed accurately via the trapezoidal rule with respect to p for any value of 9. Similarly, applying 
the trapezoidal rule in the 9 variable gives rise to high-order convergence of the integral (49), in 
view of the 7r-periodicity of the integrand. 

Remark 2. Our application of the trapezoidal rule for evaluation of 1 1 (9, uq, vq) requires use of 
equidistant samples in the p variable, which for most values of 9, do not correspond to any of 
the original grid nodes (uf ,Vm ). To address this issue our solver relies on the FFT/cubic-spline 
interpolation technique presented in [9, Section 3], which allows for fast and efficient evaluation of 
the required equidistant p samples. 
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Figure 6: Diffraction by a circular aperture: solution to the Neumann problem for an aperture of 
diameter 24A under point source illumination. The source, which is not visible here, is located to 
the left of the displayed area. The coloring on the plane is introduced for visual quality, and it does 
not represent any physical quantity. 

8.4 Type V Integral (Principal Value) 

An application of the polar change of variables mentioned in Section 8.3 to the principal- value 
Type V integral (46) results in the expression 



Here the function 4>1(p, 9) is defined by equation (52) and we have set R = r^(«o, ^o) — r i (u(p, 9), v(p, 9)). 
Equations (53) and (54) form the basis of our algorithm for evaluation of Type V integrals. 




(53) 



where Ip'^"[(j)i\(uo,vo,6) is given by the principal value integral 




(54) 



Since both 



and 



R V 

P 



are smooth functions of p and 9, it is useful to consider the expression 




(55) 



for the integral (54). This is a 1-dimensional principal value integral of the form 




(56) 



15 



Figure 7: Field diffracted by the circular aperture configuration depicted in Figure 6. From top 
left to bottom right, depiction of the diffracted field at observation screens located at distances of 
6A, 60A, 120A and 240A behind the punctured plane. A dark-spot (the Poisson shadow) can be 
observed at the center of the illuminated area in the third-to-left image. 



where v is a compactly supported smooth function. Our algorithm proceeds by evaluating this 
principal value integral by means of a trapezoidal rule algorithm with integration nodes centered 
symmetrically around x = — which, as shown in [26], yields spectral accuracy for smooth and 
periodic functions. In detail, letting x\ = (i + \)/M, after appropriate scaling into the interval 
[—1,1], our quadrature for the integral (56) is given by 



This expression provides spectrally accuracy as long as v is a smooth function of periodicity 2. Our 
Type-V integration algorithm is completed by trapezoidal integration in the 6 variable to produce 
the integral (53) with spectral accuracy. 

Remark 3. Application of the trapezoidal rule (57) to compute the integral (55) requires evenly 
spaced samples in the p variable, which, in addition, must also be symmetrically centered around 
p = 0. To obtain such samples our algorithm proceeds in two steps: 1) It uses the one-dimensional 
FFT/spline interpolation method presented in [9, Section 3] to produce evenly spaced samples of 
the integrand in the p variable, and 2) It applies an FFT-based shift (see Remark 4) to produce 
interpolated samples centered around p = 0. In view of the periodicity and smoothness of the 
function v, this procedure is highly accurate, and it is, in fact significantly faster and less memory 
intensive than the full two-dimensional spline-table construction presented in [26] — since it only 
requires storage of one-dimensional tables. 

Remark 4. Given point values v(xi) of a smooth and periodic function v on an equispaced grid 
Xi = xq + ih, samples of v on a new shifted grid x* = X{ + 5 can be obtained efficiently and with 
spectral accuracy through use of FFTs. The algorithm proceeds as follows: 1) Evaluation of the 
FFT of the data set v{x{) to produce Fourier coefficients of v(x), 2) Multiplication of each Fourier 
coefficient by an appropriate exponential, to produce the Fourier coefficients of the shifted function 
v(x + 6), and 3) Evaluation of the inverse FFT of the coefficients produced per point 2). 




(57) 
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Figure 8: Simulation of Young's experiment: diffraction by two circular apertures in a sound-hard 
plane (Neumann boundary conditions); the apertures are 24 wavelengths in diameter. The coloring 
on the plane, which is introduced for visual quality, does not represent any physical quantity. 

9 High-order evaluation of edge-patch operators 

In this section we describe our algorithms for evaluation of the edge-patch operators, namely the 
integral operators of type III, IV and VI and the differentiation operator (41). To do this, and in 
accordance with Section 4.2, we select a tensor product grid (uf 2 ,v%?) quadratically refined in v, 
which, using spatial mesh-sizes hu 2 and hf 2 in the u and t variables, is given by 

uf=uf + £hi 2 , t = \ y ...L\ 
vtf=({\ + m)hf)\ m = l...Ml (58) 

This grid is assumed to discretize a rectangle that contains H^; in view of the assumptions made 
on this set (Section 4.2) and the form of the discretization (u'' 2 ,?;™ ) we see that, while the edge 
v = is not itself sampled by this discretization, a parallel line to it, at a distance of (hf 2 /2) 2 in 
(u, v) space, is. 

9.1 Type III Integral (Regular) 

For any smooth function g defined over the interval [0, 1] which vanishes identically with all its 
derivatives at x = 1, the function g(t 2 ) can clearly be extended as a smooth and periodic function 
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Figure 9: Field diffracted by the two-hole configuration depicted in Figure 8. From left to right, 
depiction of the diffracted field at observation screens located at distances of 72A, 576A, 1728A 
and 3456A behind the punctured plane. A dark-spot can be viewed again at the center of the 
illuminated circles in the left-most image by adequately enlarging the image. 



of period 2. It follows immediately from the identity 



l fi 

+2\ 



10 V x ■ > -I ./(") 

that the trapezoidal rule approximation 

" 1 g{x) J 2 ^ , 2 . ^ 2m + 1 



dx = / g(t 2 )dt = 2 / g(t 2 )dt (59) 



r g(x) ,2^,0, 2m + 1 

/ ^~mS^(4). t m = —, l = 0,...,M-l (60) 
v m=l 



gives rise to super-algebraic convergence. Since the patch discretization («J' 2 ,Um 2 ) can be expressed 

in the form = (im 2 ) 2 , where = (i + m)hf 2 , a two-dimensional trapezoidal rule using this 
1 ;„ n i m u„i„ n ii„ 



mesh in the set H\ is super-algebraically convergent 



9.2 Partial Derivatives 

In view of the smoothness and periodicity of the function 02 (u, t 2 ), a two-dimensional interpolation 
scheme based on use of FFTs along the u variable and FCTs (Fast Cosine Transform) along the 
t variable yields spectrally convergent approximations of the function 4>2(u,t 2 ) and its derivatives. 
Our algorithm thus evaluates the derivatives required in equation (41) by performing a direct term 
by term differentiation of the resulting Fourier representations together with the expression 

dcj)2(u,t 2 ) 1 d r 2 , 

Remark 5. In view of the presence of the t = y/v denominator on the right-hand side of equa- 
tion (61), evaluation of partial derivatives of the function (f>2(u,v) with respect to v on the basis 
of a term by term differentiation of cosine expansion of the function (j)2(u,t 2 ), while yielding spec- 
trally accurate results, is less accurate near the edge than away from the edge — in close analogy 
with the well-known relative loss of accuracy around end points in Chebyshev-based numerical dif- 
ferentiation. This is why our algorithm was designed to evaluate the composite operator N W S W by 
first producing the combination "TLS^ via the rules derived for as explained in Section 7, thus 
avoiding numerical differentiation. Unfortunately, in the evaluation of the operator N w (which 
is necessary e.g. for solution of equation equation (15)), direct computation of derivatives and 
associated accuracy loss does not seem to be avoidable. 
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Figure 10: Polar changes of variables around a point close to the edge: quadratic sampling in the 
v variable, requiring off-grid interpolations for grazing angles. 

9.3 Type IV Integral (Singular) 

As in Section 8.3, we utilize a polar change of variables to resolve the Green's function singularity 
in the canonical type defined in equation (33), thus obtaining the expression 

Zf' sins [MM)= r( r H(uo,v ,p,e)— dp . )d0, (62) 
" Jo \J-^g y/v o + psm0J 

where the integrand H (no, vq, p, 9) = <f>2(uo + pcos9,vo + psm.6)^ is a smooth function of p and 
9, which vanishes for p larger than a certain constant pq. While the square-root singularity in the 
inner-integral can clearly be resolved to high-order by applying an appropriate quadratic change of 
variable, the outer integrand in 9 is not a uniformly smooth function: as detailed in Appendix B, 
it develops a boundary-layer as vo approaches 0. The analysis presented in Appendix B suggests 
a simple and efficient method for high-order resolution of this boundary layer — thus leading to 
accurate evaluation of the integral (62). This methodology, which is an integral part of our solver, 
is described in what follows. 

The aforementioned boundary-layer integration method is based on use of the change of variables 
t = \J vq + psin6*. With this change of variables equation (62) becomes 

Xl' sing [h}(uo,v ) = r I q p2 [M(^v o ,0)d9, (63) 
J o 

roo j.2 

I q p M{u ,v ,e)= H(u ,v ,—-^,9)dt. (64) 

JO SHIP 

The integral (64) is evaluated with high-order accuracy by means of a trapezoidal rule in the t 
variable, for any < 9 < n. In order to capture the boundary-layer in the outer-integral in (63), 



19 



our algorithm relies on an additional changes of variables 9 = a 2 and 9 = tt — a 2 , which lead to the 
expression 



J I p,2i ( l > 2\(uo,vo,0)de = J (l q p 2 {4> 2 ]{uQ, v , a 2 ) - I^lfe] (uo,vq, it - a 2 )) ada. (65) 

In view of the analysis presented in Appendix B, the boundary layer is confined to the interval 
[0,a*(i>o)], where a*(«o) = (^) 5 ) and we therefore decompose the a-integral in the form 

I . . .da = . . .da + J ... da. (66) 

Jo Jo Ja*(v Q ) 

For a given error tolerance, our algorithm proceeds by applying Chebyshev integration rules to 
both integrals in (66), using for the second integral a number of integration points that does not 
depend on vq, and using for the first integral a number of integration points that grows slowly 
as vq tends to zero. In practice, we have found that a mild logarithmic growth in the number 
of integration points suffices to give consistently accurate results. In view of such slow required 
growth, and for the sake of simplicity, the number of integration points used for evaluation of the 
first integral in (66) was taken to be independent of vq and sufficiently large to meet prescribed 
error tolerances; we estimate that a minimal additional computing time results from this practice 
in all of the examples considered in this paper. 

Remark 6. In order to apply the trapezoidal rule (60) for evaluation the integral of (64) we dis- 
tinguish two illustrated in Figure 10. For j < 9 < ^ we use the sampling in t provided 
by intersections with the original grid underlying V\: the 1-dimensional cubic-spline interpolation 
method introduced in section 8.3 can be used to efficiently interpolate the function JJ( f ^" , 9) at 
the needed integration points. For < 9 < | and < 9 < it, on the other hand, the i-sampling 
provided by the intersections with the original grid is too coarse. In this case, we resort to a full 
two-dimensional interpolation of the density <p2 (see Remark 7) to interpolate to a mesh in the t 
variable which, away from t = has roughly the same sampling density as that in the overall patch 
discretization. In practice a fixed number of discretization points is used to discretize all of the t 
integrals considered in the present remark. 

Remark 7. The two-dimensional interpolation method for smooth functions (f>2(u, v), which is men- 
tioned in Remark 6, proceeds by first performing a two-dimensional Fourier expansion of the func- 
tion 02 (u, t 2 ), by means of FFTs along the u variable and FCTs along the t = ^/v variable, followed 
zero-padding by a factor P (in practice we use P = 6). This procedure results in a spectral approx- 
imation of (f>2 (and, by term-by-term differentiation, of its derivatives as well) on a highly resolved 
two-dimensional grid. The final interpolation scheme is obtained by building bi-cubic spline in- 
terpolations based on function values and derivatives on each square of the refined grid. [22, p. 
195]. 

9.4 Type VI Integral (Principal Value) 

Our algorithm evaluates the principal- value edge-patch Type VI canonical integral in a manner 
similar to that used for Type IV treated in Section 9.3: introducing the local polar change of 
variables around the point r we obtain 

= /* u r HT{ua ^"- 6) , a . \ «, ( 6 7) 

•A) \ J=^% p y/vo + psmd 

\ sin / 



20 



Figure 11: Multiple scattering examples: Neumann problem on two parallel discs 24A in diameter, 
illuminated at a 45 degree angle. The coloring on the discs represents the values of the surface 
unknown 



where 

I o\ ^ R T 

H T (u o ,v o ,p,0) = 4> 2 (u + pcos8,v + psm9)—. (68) 

p 

is again, a smooth function of p and 9 which vanishes identically for p > pq. Resorting to the 
quadratic change of variables t = y/vo + p sin 8 we obtain 



Zf'HM^o)= / I*?[<h](uo,vo,9)dB. (69) 
Jo 

where the radial integral 



f°° t 2 — t 2 dt 

I q p f 1 [H(u ,v ,e)=p.v. H T (u ,v ,—?,8)-^-2, t = ^r (70) 

JO ^m V t Lq 



can be expressed in the form 



f°° v(t 2 ) 

Ip'2" = P- v - / ~2 2°^ ' w here v(t) is smooth for t > and vanishes for i large enough. 

Jo t ~ 

(71) 



'0 

Simple algebra then yields 



00 







11, f°° v(f 
\dt = p.v. / — 

t — to t + tQ \ J_ OQ t — to 



I q p f = P-v. / v(t 2 ) \-±---±-\ d t = p.v. / p^dt; (72) 
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clearly, the right hand side integral in equation (72) can be evaluated with high-order accuracy by 
means of the trapezoidal rule (57). 

Using this algorithm for evaluation of the integral (70) for any fixed value of 9 our algorithm 
for evaluation of the 9 integral (69), and thus (67), is completed, as in Section 9.3, by relying on 
(a) The quadratic change of variable 9 = a 2 and (b) The boundary-layer split (66). 

10 Parameter Selection 

A number of parameters are implicit in the algorithm laid out in Section 4, including parameters 
that relate to the overall surface patching and discretization strategies described in Section 4.2 
as well as parameters that arise in the polar integration rules introduced in Sections 8.3, 8.4, 9.3 
and 9.4. Clearly, the values of such parameters have an impact on both, accuracy for a given 
discretization density, as well as computing time for a given accuracy tolerance. A degree of 
experimentation is necessary to produce an adequate selection of such parameters for a given 
problem. Without entering a full description of the choices inherent in our own implementations, 
in what follows we provide an indication of the strategies we have used to select two types of 
parameters, namely, (a) The width of the edge patches (see Figure 1), and (b) The number of 
discretization points used in both the radial and angular directions for each polar integration 
problem for which the corresponding floating partition of unity does not vanish at the open edge, 
as illustrated in Figure 10. Similar (but simpler) considerations apply to other parameters, such as 
width of floating partitions of unity, extents of overlap between patches, etc. 

With respect to point (a) above we note that, for scattering solutions to be obtained with a fixed 
accuracy tolerance, the discretization densities must be increased as frequencies are increased, and, 
thus, the width of the edge patches can be decreased accordingly — in such a way that the number 
of discretization points in the v direction for each one of the edge patches is kept constant. This 
strategy is crucial for efficiency, since the edge patches require use of the two-dimensional interpo- 
lation method mentioned in Remark 7, which is significantly more costly than the corresponding 
one-dimensional interpolation method used in the interior patches. Use of constant number of 
f-discretization points within shrinking edge patches for increasing frequencies thus enables fixed- 
accuracy evaluation of edge-patch integrals with an overall computing cost that is not dominated 
by the edge-patch two-dimensional interpolation procedure. 

Concerning point (b) above, in turn, as mentioned in Remark 6, we make use of a fixed number 
of equispaced integration points in the scaled radial variable t (see equation (63)) for all values of the 
angular variable 9. In practice, we select the number of t-integration points to equal the maximum 
value Nt of the numbers N u and N v of points in the u-v discretization mesh that are contained in 
the 9 = and 6 = tt/2 lines, respectively, and which lie within the support of the corresponding 
floating POU. In order to preserve the wavelength sampling in the angular integral (63), finally, the 
two integrals on the right-hand-side of equation (66) are evaluated on the basis of the Clenshaw- 
Curtis quadrature rule [22] using an a discretization mesh containing points — since the length 
of half a circumference equals | times its diameter. 

11 Numerical Results 

In this section, we present results obtained by means of a C++ implementation of the algorithm 
outlined in Section 4.1, incorporating the canonical operator decompositions introduced in Sec- 
tions 5 through 7 for the operators S w , N w and N^S^, together with the high-order integration 
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N 


Dir(S w ) 


Dir(N w S w ) 


Neu(N w ) 


Neu(N w S w ) 


16 x 16 + 2 x 24 x 16 


2.4 x 10" 4 


2.5 x 10~ 4 


5.0 x 10" 4 


2.6 x 10~ 4 


32 x 32 + 2 x 48 x 32 


4.8 x 1(T 6 


4.8 x 1(T 6 


5.3 x 1(T 6 


5.2 x 1(T 6 


64 x 64 + 2 x 98 x 64 


4.7 x 10~ 8 


9.7 x 10~ 8 


4.9 x 10~ 8 


5.1 x 10~ 8 



Table 1: Scattering by a disc of diameter 3A, similar to the corresponding 24A simulation depicted 
in Figure 2. Maximum errors in the acoustic field on the square projection plate shown in the 
figure. This table demonstrates spectral convergence for all the formulations considered: doubling 
the discretization density results in orders-of-magnitude decreases in the numerical error. (The 
notation Q\ x m\ x n\ + Q2 x 771,2 x ^2 indicates that a number Q\ of patches containing m\ x n\ 
discretization points together with a number Q2 of patches containing 777-2 x "-2 discretization points 
were used for the corresponding numerical solution.) 



rules put forth in Sections 8 and 9 and the iterative linear algebra solver GMRES. Errors reported 
were evaluated through comparisons with highly-resolved numerical solutions. Computation times 
correspond to single-processor runs (on a 2.67GHz Intel core), without use of the acceleration meth- 
ods or parallelization. As mentioned in Section 4.4, application of the acceleration method [9] in 
the present context does not present difficulties; such extension will be considered in forthcoming 
work. 



Dir(SJ Dir(N w S w ) 



Disc Size 


Unknowns 


It. 


Time 


e r 




It. 


Time 


e r 




3A 


4096 


6 


58s 


1.0 x 10- 


-4 


6 


5m20s 


1.4 x 10- 


-4 


6A 


10240 


9 


3m5s 


8.2 x 10" 


-5 


6 


llml4s 


5.4 x 10" 


-5 


12A 


28672 


13 


15m21s 


1.2 x 10" 


-4 


7 


36m31s 


3.7 x 10" 


-4 


24A 


90112 


18 


2h30m 


2.8 x 10- 


-4 


7 


3h41m 


4.1 x 10' 


-4 








Neu(N w ) 






Neu(N w S w ) 




Disc Size 


Unknowns 


It. 


Time 


e r 




It. 


Time 


e r 




3A 


4096 


16 


9m21s 


1.3 x 10- 


-4 


6 


5m50s 


1.3 x 10" 


-4 


6A 


10240 


28 


36m31s 


2. x lO - 


4 


6 


llm36s 


5.8 x 10- 


-5 


12A 


28672 


49 


2h43m 


1.7 x 10- 


-4 


7 


37m04s 


5.2 x 10" 


-4 


24A 


90112 


80 


21hl0m 


1.9 x 10" 


-4 


7 


3h51m 


6.1 x 10" 


-4 



Table 2: Iteration numbers and computing times for the problem of scattering by a disc at nor- 
mal incidence. Top: Dirichlet problem. Bottom: Neumann problem. In each case, use of the 
second- kind combined operator N W S W gives rise to significantly smaller iteration numbers than the 
corresponding first kind formulation. In the case of the Neumann problem, the reduction in iter- 
ation numbers results in substantially improved computing times. Note: all reported computing 
times correspond to non- accelerated single-processor runs. Dramatic reductions in computing times 
would result from use of the acceleration method [9] — see e.g. the recent contribution [8] for the 
closed-surface case. 
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11.1 Spectral convergence 

We demonstrate the spectral properties of our algorithm through an example concerning a canonical 
geometry, namely, the unit disc 

x 2 + y 2 < 1, z = 0. (73) 

For this surface we utilize three coordinate patches (see Figure 1), including a large central patch 
given by equations {x(u,v) = u, y{u,v) = v, z(u,v) = 0}, and two edge patches parametrized 
by the equations {x(u, v) = (1 — v) cos u, y(u, v) = (1 — v) sinu, z(u, v) = 0} for values of u and 
v in adequately chosen intervals. The two edge-patches overlap as illustrated in Figure 1, and 
their width is defined by the range of the v variable, which, in accordance with Section 10, is 
reduced as the frequency increases. With reference to equation (12), the integral weight is set to 
uj = \J\ — x 2 — y 2 . 




Figure 12: Poisson-spot phenomenon. Left: cross-sectional view on a of the diffraction pattern pro- 
duced by a disc 24A in diameter in three dimensional space (Dirichlet problem, normal incidence). 
Right: Diffraction by an arc of length 24A in two dimensional space (Dirichlet problem, normal 
incidence). Note that only in the three-dimensional case does a "Poisson cone" and corresponding 
"Poisson spot" develop in the shadow region. 

The sound-soft (Dirichlet) problem can be solved by means of either the first-kind equation (14), 
or the second-kind equation (16) which, in what follows, are called Dir(S w ) and Dir(N w S tJ ), re- 
spectively. The sound-hard (Neumann) problem, similarly, can be tackled by means of either the 
first-kind equation (15) or the second-kind equation (17); we call these equations Neu(N w ) and 
Neu(N w N aj ), respectively. Table 1 demonstrates the high-order convergence of the solutions pro- 
duced by our implementations for each one of these equations on a disc of diameter 3A; clearly 
errors decrease by orders of magnitude as a result of a mere doubling of the discretization density. 

11.2 Solver performance under various integral formulations 

In this section we demonstrate the performance of the open-surface solvers based on use of the 
operators Dir(S^) and Dir(N w S w ) for the Dirichlet problem, as well as the operators Neu(N w ) and 
Neu(N w S w ) for the Neumann problem. We base our demonstrations on two open surfaces: a disc 
and a spherical cavity defined by 

x 2 + y 2 + z 2 = l, z > cos(0 o ), (74) 
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Dir(S^) Dir(N^) 



Spherical Cavity Size Unknowns 


It. Time e r 


It. Time e r 


3A 9344 


17 7ml6s 1.8 x 10~ 4 


13 lhl8m 4.5 x 10~ 4 


9A 84096 


39 4h06m 2.1 x 10~ 4 


24 13h20m 2.9 x 10~ 4 


18A 336384 


65 57h48m 4.0 x 10~ 4 


43 124h 1.4 x 10" 4 


Neu(Nj Neu(N w S w ) 


Spherical Cavity Size Unknowns 


It. Time e r 


It. Time e r 


3A 9344 


57 lh24m 1.2 x 10~ 4 


13 lh20m 8.8 x 10" 4 


9A 84096 


243 52hl4m 2. x 10" 4 


24 13h21m 5.6 x 10" 4 


18A 336384 


> 600 


43 124h 3.1 x 10" 4 



Table 3: Iteration numbers and computing times for the problem of scattering by the spherical 
cavity defined by equation (74) and depicted in Figures 4 and 5. Top: Dirichlet problem. Bottom: 
Neumann problem. Reductions in numbers of iterations and computing times occur as detailed in 
the caption of Table 2, but, owing to the rich multiple scattering phenomena that arise within the 
cavity, the iteration numbers are significantly higher, in all cavity cases, than those required for 
the corresponding disc problems. 



where 9$ denotes the cavity aperture. For the examples discussed here we set #o = if> an< ^ we 
made use of the weight function u = yj z — zq where zq = cos(#o)- 

For both geometries, computational times and accuracies at increasingly large frequencies are 
reported in Tables 2 and 3. In all tables the acronym It. denotes the number of iterations required 
to achieve a relative error (in a screen placed at some distance from the diffracting surface) equal 
to "e r " (relative the the maximum field value on the screen), and "Time" denotes the total time 
required by the solver to evaluate the solution. As can be seen from these tables, the equation 
Neu(N w ) requires very large number of iterations for the higher frequencies. The computing times 
required by the low- iteration equation Neu(N w S w ) are thus significantly lower than those required 
by Neu(N w ). The situation is reversed for the Dirichlet problem: although the equation Dir(N w S tJ ) 
requires fewer iterations than Dir(S w ), the total computational cost of the low-iteration equation is 
significantly higher in this case — since the application of the operator in Dir(S w ), which, fortunately, 
suffices for the solution of the Dirichlet problem, is substantially less expensive than the application 
of the operator in Dir(N w S w ). As it happens, at high-frequency, the bulk of the computational 
time used by our solver is spent on interior-patch work: application of the acceleration method 
of [9] (see also [8]) would therefore reduce dramatically overall computing times for high-frequency 
problems. 

Figures 2 and 3 display three dimensional renderings of patterns of diffraction by the disc 
under normal incidence with Neumann boundary condition, and under horizontal incidence with 
Dirichlet boundary condition. Corresponding images for the spherical-cavity problem are presented 
in Figures 4 and 5; note the interesting patterns of multiple-scattering and caustics that arise in 
the cavity interior. 



25 



Figure 13: Dirichlet problem on an array of 8 x 8 discs of diameter 6A. (Overall diameter: 96. 6A; 
192 patches used.) The coloring on the discs represent the values of the surface unknown ip. 

11.3 Miscellaneous examples 

This section presents a variety of results produced by the open-surface solver introduced in this 
paper, including demonstration of well known effects such as the Poisson spot, and applications in 
classical contexts such as that provided by the Young experiment. 

11.3.1 Poisson spot 

As mentioned in the Introduction, the experimental observation of a bright area in the shadow of 
the disc, the famous Poisson spot, provided one of the earliest confirmations of the wave-theory 
models of light. The Poisson spot is clearly visible in the diffraction patterns presented in Figure 2 
and the left portion of Figure 12. The left portion of Figure 12 displays a slice of the total field 
around the disc along the x — z plane, which gives a better view of the Poisson-spot phenomenon: 
the "Poisson cone" is clearly visible in this figure. Interestingly, this phenomenon does not occur 
in the two-dimensional case. This is demonstrated in the image presented on the right portion 
of Figure 12: the two-dimensional diffraction pattern arising from the flat unit strip (which was 
obtained by the solver presented in [10]) gives rise to a dark dark shadow area which does not 
contain a diffraction spot. 

11.4 Babinet's principle, apertures and Young's experiment 

For a flat open surface T contained in a plane II one may consider the corresponding problem of 
diffraction by the complement r c = II \ T of T within II. As is well known, the diffraction pattern 
resulting from r c can be computed easily, by means of the Babinet principle (see Appendix C 
and, in particular, equation (93)), from a corresponding diffraction pattern associated with the 
surface T. (For ease of reference, a derivation of the Babinet principle for scalar waves is presented 
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Figure 14: Neumann problem for an array of 8 x 8 circular apertures of diameter 6A. The diffracted 
field depicted in this figure was produced by means of Babinet's principle from the diffraction 
pattern displayed in Figure 13. As in Figure 6, the coloring on the plane is introduced for visual 
quality, and does not represent any physical quantity. 

in Appendix C.) In what follows we present three applications of the Babinet principle, namely, 
the diffraction by a circular aperture, the Young phenomenon, and diffraction across an array of 
apertures (in Section 11.4.1). 

As our first application of Babinet's principle, in Figure 6 we present the field diffracted by 
a circular aperture which is 24A in diameter. The incident field for this image was taken to be 
a point source located at the point (0,0,-10), outside the region displayed on the figure. In 
Figure 7, we display the total field on screens located behind the aperture at varying distance from 
the punctured plane. Interestingly, under some configurations a dark spot appears in the center 
of the bright area, in full accordance with Arago's prediction that a 'Poisson shadow' must exist. 
As our second application of Babinet's principle, in Figure 8 we present the field diffracted by a 
pair of nearby circular holes in an otherwise perfectly sound-hard plane, under normal plane-wave 
incidence: this is a setup of the classical Young experiment. This diffraction pattern was produced, 
by means of Babinet's principle, from a corresponding solution of a Dirichlet problem for two 
coplanar discs in space. As in Young's experiment, interference fringes arise: these can be seen 
clearly on the right-most image in Figure 9. Again, sharp dark spots at the center of the circular 
illuminated areas can be seen in the leftmost image in Figure 9. 

11.4.1 Arrays of scatterers/apertures 

Geometries consisting of a number of disjoint open-surface scattering bodies can be treated easily 
by the solver introduced in this paper — since the decomposition in patches inherent in equation (19) 
is not restricted to sets of patches representing a connected surface. The solution for the two-disc 
diffraction problem presented in the previous section, for example, was obtained in this manner. 
In what follows we provide a few additional test cases involving composites of open surfaces. 
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Figure 1 1 presents the solution of a problem of scattering by two parallel discs illuminated at an 
angle of j, with Neumann boundary conditions. The beam reflected by the bottom disc, which can 
clearly be traced onto the upper disc, gives rise to a bright area in the projection screen behind that 
disc. Two final examples concern an array of 64 discs and the corresponding array of 64 circular 
apertures on a plane — where each disc is 6A in diameter and the discs are separated by 3A spacings, 
for a total array diameter of 96. 6A. The corresponding diffracted fields are presented in Figures 13 
and 14. The solution of the Dirichlet problem for the 64 disc array was obtained in 23 GMRES 
iterations on a 192 patch geometry representation. 

12 Conclusions 

We have introduced a new set of integral equations and associated high-order numerical algorithms 
for the solution of scalar problems of diffraction by open surfaces. The new open-surface solvers 
are the first ones in the literature that produce high-order solutions in reduced number of GM- 
RES iterations for general (smooth) open surfaces and arbitrary frequencies. The second-kind 
formulation Neu(N w S tJ ) is highly beneficial in the context of the Neumann problem, as it requires 
computing times that are orders-of-magnitude shorter than those required by the alternative hy- 
persingular formulation N(N W ). Such gains do not occur for the Dirichlet problem: the proposed 
solver produces high-order solutions to Dir(S^) in very short computational times, and the gains in 
iteration numbers that result from use of the formulation Dir(N w S w ) do not suffice to compensate 
for the significantly higher cost required for evaluation of the operator N w . With appropriate par- 
allelization and acceleration, fast and accurate solutions should be achievable by these algorithms 
for very large structures, thus providing a robust numerical workbench for the solution of classical 
diffraction problems by open surfaces. 

Acknowledgments. The authors gratefully acknowledge support from AFOSR, NSF, JPL and 
the Betty and Gordon Moore Foundation. 

A Expression of the operator N w in terms of tangential derivatives 

In this section, we provide a proof of Lemma 1. 

Proof. It suffices to show that the operators on the left- and right-hand sides of equation (35) 
coincide when applied to any smooth function ip defined on T. Following the derivation in [12] for 
the closed-surface case, we define the weighted double layer operator 

D w M(r) = jf ^^(r'Mr')dS', (75) 

and we evaluate the limit of its gradient as r tends to T. For r outside T, (75) can be expressed in 
the form 

D w M(r) = -div ^G k (r,r')4,(r')Lu(r')n r ,dS'. (76) 

Thus, using the identity 

curl curLA = -AA + Vdiv A, 
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we obtain 



VD^](r) = fc 2 f G k (r,r'W(r')u(r')n rl dS' 

JF , (77) 
-curl curl J G k (r, r')^(r')n r ,dS' . 

But, for r outside T we have 

curl ! G k (r,r')^(r')n r >dS' = 

Jr (78) 
[n r ^^r>(r')V*,G fc (r,r')] dS' , 

where [ • , • ] and V s denote the vector product and surface gradient operator, respectively. Inte- 
grating by parts the surface gradient (see e.g. [12, eq. (2.2)]) and noting that the boundary terms 
vanish in view of the presence of the weight u>, we obtain 



curl J G k (r,r')^(r')n r ,dS' = 
-j G k (ry) [rv,V s (^(r>(r'))] dS'. 
In the limit as r tends to an interior point in T we therefore obtain the expression 

VD w M(r) = k 2 ^G k (r,r'W(r')u(r')n r ,dS' 
+p.v.^ [V r G fe (r,r'), [n r ,, V s (V(r>(r'))]] dS' (r G T), 



(79) 



(80) 



in terms of a principal value integral, for the surface values of the gradient of the double layer 
operator D w . Taking the scalar product with n r on both sides of (80) now yields the desired result: 
equation (35). 

□ 

B Boundary-layer character of the inner integral in equation (62) 

In order to demonstrate the difficulties inherent in the numerical evaluation of the outer integral 
in equation (62) we consider the integration problem 



I p (v ,p ,6)dd, (81) 

io 

in which the (uo, ^-dependent inner integral in (62) is substituted by the vo-dependent integral 

f°° ~ do 
I p (v ,p ,e)= Hp() (p,9)-=IL= , 0<e<7r. (82) 

J__2L V vo + p sin 

sin 6 

Here 

\ri>2. (83) 
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so that, in the present example, po is the polar-integration radius. (The inner integral in (62) varies 
smoothly with uq, and, thus, the uq dependence does not need to be built into the present analogy.) 
The integral I p (vq, po, 9) is given by 



I P (vo,po,9) 



7 = . I 1 * if sin6><^ 

V'-'O+PO sinS+vvo— po sin6 — po 



Vsin 6 Po 



(84) 



Clearly, as vo tends to zero, I(vq, po,9) becomes increasingly singular (as demonstrated in the left 
portion of Figure 15): in view of the last equation on the right hand side of (84), we have 

]hnJ p (v o ,p o ,0)=(—LJ). (85) 



"o->-o \vsin 

To treat the singularity in (85) we introduce quadratic changes of variables in the 9 integration in 
equation (81) — of the form = a 2 in the interval [0, and 9 = II — a 2 in the interval [f , vr]). As a 
result of these operations we obtain bounded integrands: for example, the integrand resulting from 
the first of these changes of variables is J(vo, po, a) = aI(vo, po, a 2 ), which is a bounded function 
of a. This integrand is depicted on the right portion of Figure 15; clearly J(vo, po, a) develops a 
boundary layer as vo tends to zero. 

The two changes of variables mentioned above result in integrals over the domain [0, y^]), and 
in both cases boundary layers result at and around a = 0. To resolve these boundary layers we de- 
compose the integration interval into two sub-intervals, namely [0, a* (vo, po)] and [a*(vo, po), \/f-]- 
Here, for a given value of vq, the point a*(vo,po) is chosen to lie to right of the coordinate for 
which the peak occurs in the right portion of Figure 15, in such a way that the slope of the function 
J( v 0i Po, a ) as a function of a at a = a*(vo,po) remains constant as vo approaches zero — with 
a slope that equals a certain user-prescribed constant value. In practice we have found that the 
integral to the right of the point a = a*(vo,po) can be performed, with fixed accuracy, by means 
of a number of discretization points that grows very slowly as vq —> 0. (In our implementations we 
typically use a number of discretization points to evaluate this integral that remains constant for all 
required small values of vo-) The evaluation of the integral on the left of the point a = a*(vo,po) 
with fixed accuracy requires a number of discretization points that does grow somewhat faster, 
as vo — > 0, than the one on the right, but, we have found in practice that the latter integral can 
be obtained with fixed accuracy by means of a number of discretization points that grows only 
logarithmically with vq as vo — > 0. 

To obtain an approximate expression for a*(vo,po) we note that, since 



it ^ V v o + Po sin a 2 I . ,v 

J (vo, po,a) = a ^ tor a>,/arcsm( — ), (86) 

since 2 y p 

for a 1 and for sufficiently small values of vq (such that the inequality constraint in (86) is 
satisfied) we have a 2 ~ sin a 2 , and thus letting rj = ^ , 

J(v ,p ,a) ~ </pofn(a) where f v (a) = " • (87) 

It follows that, for a given constant C > 0, the fixed-slope point for which /' (a^) = —C is 
approximately given as a root of the equation 

77 OL^ eft 1 

i 1 = -C, or equivalents, ^ + ^ - — = 0. (88) 

a 2 ^J V + a 2 V V 2 C 2 
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Clearly, thus, an approximation of the quantity a 2 can be obtained, in closed form, as a root of 
a certain polynomial of degree three. A Taylor expansion of the resulting root as a function of r/ 
around rj = shows that 



a,. 



as i] — > 0, 



(89) 



and, since rj 

by 



it follows that the constant slope point a*(vo,po) is given, for each constant po, 



a*(v ,po) ~ u 3 as v 



0. 



(90) 



Since the function .ff = H(uq, Vq, p, 9) in equation (62) is modulated by a smooth windowing 
function that is akin to the "discontinuous window function" H po in equation (84), it is reasonable 
to expect that the inner p-integral in (62) gives rise to a a-integrand which develops a similarly 
behaved, albeit smoother, boundary-layer. We illustrate this in Figure 16 (the left portion of which 
should be compared to the right portion of Figure 15), which displays the function 



J(v ,p ,a) = a 



H(v ,p, a 2 



POP 



\/vq + 



po sin or 



where 



P 



H po (p,9) = W(^), W(p) 
Po 



e 1 -p 2 , p < 1 
0, p>l. 



(91) 



(92) 



It follows that the ^-integration strategy outlined above in this section for the function / (based on 
the change of variables 9 = a 2 and partitioning of integration intervals at the point a = a*) applies 
to the integrand given by the inner integral in equation (62): this strategy is incorporated as part 
of our algorithm, and is thus demonstrated in the numerical examples presented in Section 11. 



A 



v0=0.1 
-v0=0.05 
-v0=0.01 



0.5 



1.5 



1.5 



0.5 



N 



0.2 



0.4 



0.6 



0.8 



- v0=0.1 
— v0=0.05 
— v0=0.01 



1.2 



1.4 



Figure 15: Left: boundary layer for I(vq, po, 9) on the interval [0, |] 
J{vq, po, a) = aI(vo, po, a 2 ), where a G [0, 



Right: quadratic regularization 



C Babinet Principle for Acoustic Problems 

As mentioned in Section 11.4, for ease of reference, in this appendix we present a derivation of the 
Babinet principle for scalar waves; see also [7]. Let T be an open (bounded) flat screen which lies 
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Figure 16: Left: numerical values of J{vq, pq,ol) for the smooth function H po (p,8) given in equa- 
tion (92). Right: normalized view on the interval [0, a*(vo)] for various values of vq. 

in the z = plane, and let u l and up denote the incident wave (with sources contained in the semi- 
space z < 0), and the corresponding field scattered by T under the Dirichlet boundary condition 
given by equation (1) with / = — u?. The corresponding total field is denoted as u-f = u l + up. 

Calling T c the complement of T in the z = plane, let vp c denote the field scattered by T c 
under Neumann boundary conditions, and let i>p c denote the corresponding total field. We now 
establish the Babinet principle that relates the Dirichlet-screen problem to the Neumann aperture 
problem, namely 

+ vfc = u l for z > 0, (93) 

with an associated formula, given below, for z < 0. The corresponding Babinet principle relating 
the Neumann-Screen problem to the Dirichlet-aperture problem follows similarly. 

In the Dirichlet-screen/Neumann-aperture problem the total field t>p c satisfies Up c = u % on T — 
since vf c = u l + up c , and since up c is an odd function of z (as it equals a double-layer potential 
with source density on T c ). Thus, defining w(x,y,z) = v^ c (x,y,z) for z > 0, and w(x,y,z) = 
Vfc(x,y, —z) for z < 0, we see that w satisfies the following properties 

— The boundary values of w on T satisfy i0| r = uj ; 

— w is a radiative solution in M 3 (since Up c is radiating behind the screen) 

- w is continuous across T c (by definition) and its normal derivative across T c is continuous 
(since t>p c satisfies homogeneous Neumann boundary conditions on T c ). 

It follows that w = —up (by uniqueness of solution to the Dirichlet problem on T) and, therefore, 
that vf c = —up for z > or, in other words, that vp c = —up — u l for z > — and thus equation (93) 
follows. Since, as stated above, up c is an odd function, its values in the region z < follow by 
symmetry. 
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